Ground-State Structures of Ice at High-Pressures 

Jeffrey M. McMahon^Q 

^Department of Physics, University of Illinois at Urbana- Champaign, Illinois 61801, USA 

(Dated: June 13, 2011) 

Abstract 

Ah initio random structure searching based on density functional theory is used to determine the 
ground-state structures of ice at high pressures. Including estimates of lattice zero-point energies, 
ice is found to adopt three novel crystal phases. The underlying sub-lattice of O atoms remains 
similar among them, and the transitions can be characterized by reorganizations of the hydrogen 
bonds. The symmetric hydrogen bonds of ice X and Pbcm are initially lost as ice transforms to 
structures with symmetries Pmc2i (800 - 950 GPa) and P2i (1.17 TPa), but they are eventually 
regained at 5.62 TPa in a layered structure C2/m. The P2i — t- C2/m transformation also marks 
the insulator-to-metal transition in ice, which occurs at a significantly higher pressure than recently 
predicted. 

PACS numbers: 64.70.K-, 62.50.-p, 71.30.+h, 96.15.Nd 
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The behavior of H2O at high pressures is of fundamental important for both condensed 
matter and planetary physics [H |2]. This can be attributed to its substantial abundance in 
the universe, and the fact that a significant fraction of it exists in ice form at high pressures 
in planetary interiors. In our solar system alone, for example, Uranus and Neptune consist 
largely of H2O, ammonia, and methane ice mixtures up to 800 GPa, and the cores of Saturn 
and Jupiter likely contain ice components at pressures of approximately 800 GPa - 1.8 TPa 
and 4-5 TPa, respectively [3]. Despite this importance, very little is known about the 
behavior of solid H2O (ice) under such extreme conditions. The primary reasons for this are 
that experiments have thus far only reached 210 GPa [4J and, until recently, theoretical and 
computational methods have not existed to reliably predict crystal structures with little or 
no a priori information. 

The known phase diagram of H2O is extremely rich. 15 thermodynamically stable phases 
of ice have been observed experimentally [5], ice XV only recently [6]. The highest pressure 
phase experimentally observed is ice X, which is obtained from a phase-transition from ice 
VII (or VIII, depending on the temperature) near 44 GPa [7j. In this phase, the O atoms 
form a body-centered cubic sub-lattice and the H atoms adopt symmetric positions between 
them at pressures near 110 - 120 GPa [8J, and so this phase is also referred to as symmetric 
ice. Because of this the distinction between covalent bonds and hydrogen bonds is lost, as 
is thus the molecular form of H2O, resulting in an atomic solid. Recent lattice dynamical 
calculations using density functional theory (DFT) suggest that the symmetric ordered form 
of ice X is only stable from 120 - 400 GPa 0. Near 300 - 400 GPa a lattice instability 
occurs, resulting in a transition to a crystal phase with Pbcm symmetry (Hermann-Mauguin 
space-group symbol) [9], as first predicted by Benoit et al. in 1996 via constant pressure 
molecular dynamics simulations fTO^. This phase results from a compromise between the 
packing efficiency of the O atoms and a preservation of the symmetric hydrogen bonds, 
resulting in a distorted hexagonal close-packed (hep) sub-lattice of O atoms. At even higher 
pressures, which is of fundamental importance to planetary physics, for example, very little 
is known. However, a number of intriguing possibilities have been proposed, perhaps the 
most interesting being an insulator-to-metal transition pT| [T2] . 

In this Letter, the recently proposed ah initio random structure searching (AIRSS) 
method of Pickard and Needs to predict crystal structures jl3j, combined with DFT, is 
used to determine the ground-state structures of ice at high pressures. In this method, a 
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number of random configurations are each relaxed to the ground state at constant pressure 
(see below). After enough trials, a good sampling of the configuration space is obtained and 
the ground-state structure(s) can be identified. This method has been used to successfully 
predict the ground-state structures of a number of systems most recently including 
those of atomic metallic hydrogen p!5] . 

Calculations were performed using the Quantum ESPRESSO DFT code [16]. Norm- 
conserving TrouUier-Martins pseudopotentials [17] were used for all calculations. For O, 
core radii of 1.25, 1.25, and 1.4 a.u. were used for the s, p, and d components, respectively. 
For H, a core radius of 0.8 a.u. was used for the AIRSS and then decreased to 0.3 a.u. for 
recalculating detailed enthalpy vs pressure curves. For all calculations, the Perdew-Burke- 
Ernzerhof generalized-gradient approximation exchange and correlation functional [18] was 
used. See Ref. [19] for justifications for these approximations. A plane-wave basis set with 
a cutoff of 120 Ry was used for the AIRSS and then increased to 350 Ry for recalculating 
enthalpy curves. For Brillouin-zone sampling, 8^ k-points were used for all calculations, 
except for Cmcm, Cmca, P4:2/nnm, and P2i/m (see below) for which 12'^ were used to 
recalculate enthalpy curves. The high-accuracy cutoff and k-point sampling were found to 
give a total convergence in energy to better than 1 mRy/H20 and the total pressure to 
approximately 1 GPa for each structure. Phonons were calculated using density functional 
perturbation theory as implemented within Quantum ESPRESSO, and were converged to a 
similar level of accuracy as the DFT calculations. 

Random structures were constructed by generating random unit-cell translation vectors, 
renormalizing the volume, and choosing random H2O configurations (positions and orienta- 
tions). Constant pressure geometry relaxations were performed at 0.5, 1, 1.5, 2, 3, and 5 
TPa for unit cells containing 4 H2O units, and then additional relaxations were performed 
at 1 and 2 TPa using unit cells containing 6 and 8 H2O units. (However, the latter searches 
only revealed a couple of additional structures - see below and Ref. [H].) It is important to 
realize that searches over unit cells of their factors are implicitly included in these calcula- 
tions - i.e., those with 1, 2, or 3 H2O units. While structures with unit cells containing 5, 7, 
or more II2O units are certainly possible, it is reasonable to suspect that they are unlikely 
based on comparisons with the other predicted high-pressure phases of ice, such as Pbcm 
[TU] and the recently proposed Cmcm and Pbca structures [12] (all found via simulations 
capable of generating unit cells with up to 16 II2O units). Typical relaxations included 
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FIG. 1: (color online). Enthalpies of the ground-state structures of ice relative to Pbcm, not 
including lattice zero-point energies. Note that the enthalpies vs pressure are nearly linear from 
2.5 - 5 TPa (not shown). 

up to 175 random structures at each pressure considered, which appeared to be enough to 
generate both the lowest-enthalpy structure and higher-enthalpy metastable ones multiple 
times. 

As a confirmation of the validity of this method to reliably predict the high-pressure 
phases of ice, at 500 GPa 10% of the relaxations revealed Pbcm as the lowest-enthalpy 
structure. This is consistent with previous predictions P HU], indicating that AIRSS as 
described above is indeed capable of generating (presumably correct) ground-state structures 
with a high probability and that the searches are exhaustive, as well as supports the ice X 
— > Pbcm transformation. 

After performing the AIRSS, each structure within approximately 25 mRy/H20 of the 
lowest-enthalpy one found was considered for further investigation. Detailed enthalpy vs 
pressure curves were calculated for these structures by performing additional constant- 
pressure geometry optimizations while keeping symmetries fixed; Fig. [T} 

Near 800 GPa, Pbcm becomes unstable relative to two additional structures, Pmc2i 
and Pbca. Such instability is expected, as it has been demonstrated that Pbcm develops a 
dynamic instability near 760 GPa in the (1/2, 0, 0) phonon mode [12]. Pbcm and Pmc2i 
are both shown in Fig. [2| and Pbca is shown in Refs. [121 HH] . Pmc2i and Pbca are both 
similar to Pbcm. For example, in the perspective of Fig. [2] the O atoms remain close to their 




4 




FIG. 2: (color online). Ground-state structures of ice at 1 TPa. (left) Pbcm and (right) Pmc2i. 

distorted hep sub-lattiee positions [TU]. However, the H atoms are shifted away from their 
symmetric 0-H-O positions. In Pbca, a small distortion of the H atoms occurs in alternating 
directions, leaving them close to tetrahedral sites and the hydrogen-bond network intact. In 
Pmc2i, on the other hand, a reorganization of the hydrogen-bond network occurs. In the 
perspective of Fig. [2| every O atom in Pbcm is connected to both its vertical and horizontal 
neighbors via symmetric hydrogen bonds. In Pmc2i, this bonding is only retained in every 
other column of O atoms, but without symmetric hydrogen bonds. The other O atoms 
become disconnected from their vertical neighbors, and instead become hydrogen bonded 
(and not symmetrically) with O atoms out of the plane. 

It should be noted that Pbca was recently proposed as a likely candidate for the ground- 
state structure of ice from 760 GPa to 1.25 TPa [12]. However, Fig. [T] shows that it is only 
competitively stable with Pmc2i from 800 - 925 GPa. It is certainly possible that Pbca is 
a stable phase of ice in this narrow pressure range. However, given that it is only a slight 
distortion of Pbcm and less than 1 mRy/H20 more stable, a more likely scenario is that it 
is a transitional structure that occurs during the Pbcm — ?■ Pmc2i transformation (this is 
further suggested by lattice dynamics calculations that show Pbca is unstable near 800 GPa 
- see Ref. [I9]). In any case, 800 GPa (approximately) marks a transition away from the 
symmetric hydrogen bonds seen in ice X and Pbcm, which are shown below to be regained 
at much higher pressures. 

Pmc2i remains stable until 1.3 TPa. It then becomes relatively unstable towards a struc- 
ture with P2i symmetry, which is shown in Fig. |3} Comparison of Figs. [2] and [s]^ left) shows 
that (in the plane of each figure) the O atoms continue to remain close to their distorted hep 
positions. Into the plane, however, P2i undergoes a noticeable compression and slight dis- 
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view. 



FIG. 4: (color online). C2/m at 5 TPa. The perspectives are similar to those shown in Fig. [sj 
Note that alternating O atoms are out of the plane relative to one another. 

tortion relative to Pbcm (not shown). Moreover, it can be seen that a further reorganization 
of the hydrogen-bond network occurs. O atoms in every other column continues to remain 
connected to their vertical neighbors, but the hydrogen bonds distort slightly outwards in 
alternating directions. The hydrogen bonds of the other O atoms, on the other hand, re- 
arrange more significantly, and each O atom becomes connected to either one neighboring 
column or the other, also in an alternating fashion. 

P2 1 remains the lowest-enthalpy structure up to the highest pressure considered in this 
work of 5 TPa. However, another competitive structure was also found in this pressure 
range, C2/m. In fact. Fig. [l] shows that the relative enthalpy difference between C2/m and 
P2i slowly decreases with increasing pressure. A linear extrapolation of the enthalpy vs 
pressure curves near 5 TPa (which should be quite accurate) indicates a transition pressure 
of approximately 5.29 TPa. As can be seen in Fig. |4| C2/m is a layered structure, where 
each layer consists of two sets of O atoms. Additionally, a slight shear deformation of the 
layers occurs, leaving the O atoms close to, but slightly displaced from their distorted hep 
positions. (Note that a less stable structure without the shear deformation, C2/c, was also 
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FIG. 5: (color online). Metastable structures of ice at 2 TPa. (left) P-Sml and (right) Fddd. 

found - see Ref. [H].) Moreover, symmetric hydrogen bonds are seen to be regained, which 
connect the O atoms within each layer. 

Comparison of the O atoms in Fig. |3]( right) with |4]^left) indicates that a transition to- 
wards the layered Cljvn structure is evident in P2i (and in fact even Pbcm, which looks 
similar but less compressed, as discussed above). Given these results, it is possible to un- 
derstand the phase transformations that ice undergoes at high pressures. In all structures, 
the O atoms remain close to their distorted hep positions and compress along a preferred 
axis (e.g., into the plane of Fig. [2]). This compression leads to structures without symmetric 
hydrogen bonds, Pmc2i and P2i, obtained via reorganizations of the hydrogen-bond net- 
work. Continued compression eventually results in a layered structure where the symmetric 
hydrogen bonds are regained. 

A number of structures with higher enthalpies were also generated during the AIRSS, as 
evident from Fig. [!} It has recently been demonstrated (at medium pressure) that such 
metastable structures of ice are both experimentally realizable [20] and predictable via 
AIRSS [21]. It is thus quite possible that those shown in Fig. [l] may actually occur in 
nature, and therefore worthwhile to very briefly mention a couple of them here. Brief de- 
scriptions of the other ones can be found in Ref. [T9j . Two of the most interesting metastable 
structures found were P-3ml and Fddd, as shown in Fig. [5j In both structures, the O atoms 
once again remain near distorted hep positions. However, their hydrogen bonds are quite 
different. In P-3ml, for example, 2/3 of the O atoms have 3 strongly-coordinated H atoms 
and the other 1/3 have none, but 6 weakly-coordinated ones instead. Fddd, on the other 
hand, forms a string-like structure though its hydrogen bonds. (Note that below 1.6 TPa, 
Fddd distorts into a lower symmetric form, Fdd2.) 

The results presented above were for static lattices. However, the light hydrogen mass 
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causes the phases of ice at high pressures to have large zero-point energies (ZPEs) that 
must be estimated in order to determine the most stable ground-state structures. ZPEs 
were neglected during the AIRSS, but their impacts were estimated afterwards using the 
harmonic approximation: -Ezpe = J du F{uj)huj/2, where F{uj) is the phonon density of 
states. F{uj) was calculated using a 2^ grid of q-points in the Brillouin zone, which is 
estimated to be sufficient to converge ZPE differences between structures to within a few 
percent. 

Reference [12] shows that the ZPEs are quite large, increasing from approximately 67 
mRy/H20 at 400 GPa to 127 mRy/H20 at 5 TPa. Despite such large values, ZPE differences 
between the structures are relatively small, in all cases within a few mRy/H20. While this 
energy scale is not enough to change the ordering of the structures, it is enough to affect 
precise transition pressures, in some cases. For example, P2i is found to have a lower ZPE 
than both Pmc2i and C2/m, causing the corresponding transition pressures to be shifted 
to 1.17 and 5.62 TPa, respectively. The ZPE difference between Pbcm and Pmc2i, on 
the other hand, is found to be practically negligible, resulting in a shift of the transition 
pressure higher by less than 50 GPa. Note that these estimates neglect the impact of zero- 
point pressure, which given the small differences in ZPE between the structures should be 
relatively minor. 

One of the most intriguing suggestions regarding high-pressure ice is its metallization 
[llj . Pursuant to this, the electronic density of states for each structure shown in Fig. [l] 
was investigated. Except for Fddd and Cmm2, all of the structures eventually become 
metallic with increasing pressure. However, Pmc2i was found to be insulating over its 
entire range of stability as was P2i up to 4.7 TPa. Furthermore, given that DFT typically 
underestimates band gaps, the actual insulator-to-metal transition in P2i likely occurs at 
a much higher pressure. Although, metallization in C2/m occurs at a pressure much lower 
than the P2i — )■ C2/m transformation (in fact, it was found to be metallic at all pressures 
considered). It can therefore concluded that this phase transformation at 5.62 TPa also 
marks the insulator-to-metal transition in ice, which is a significantly higher pressure than 
the recent prediction of 1.55 TPa [T2] . 

In conclusion, AIRSS was used to determine the ground-state and metastable structures 
of ice at high pressures. The predicted transformation sequence is ice X — )■ Pbcm (300 - 
400 GPa [IDJ) ^ Pmc2i (800 - 950 GPa) P2i (1.17 TPa) C2/m (5.62 TPa), where 
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transition pressures have been indicated in parenthesis. No additional competitive structures 
were found during the AIRSS, even at 5 TPa. Thus, any phases existing beyond C2/m do 
so at much higher pressures, which are clearly outside of the range of experimental and 
most astrophysical applications. The P2i — > C2/m transformation was demonstrated to 
mark the insulator-to-metal transition in ice, which is also beyond pressures found inside 
even many giant planets, such as Jupiter, and it can therefore be concluded that the ice 
components in them remain insulating. Along with the recent work elucidating the high- 
pressure high-temperature phase diagram of water [22], the results presented herein provide 
a comprehensive picture of the high-pressure phase diagram of H2O. 
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